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INTRODUCTION 


Solution of the wave equation using techniques such as finite difference or finite 
element methods can model elastic wave propagation in solids. This requires mapping the 
physical geometry into a computational domain whose size is governed by the size of the 
physical domain of interest and by the required resolution. This computational domain, in 
turn, dictates the computer memory requirements as well as the calculation time. Quite 
often, the physical region of interest is only a part of the whole physical body, and does not 
necessarily include all the physical boundaries. Reduction of the calculation domain 
requires positioning an artificial boundary or region where a physical boundary does not 
exist. It is important however that such a boundary, or region, will not affect the internal 
domain, i.e., it should not cause reflections that propagate back into the material. This paper 
concentrates on the issue of constructing such a boundary region. 


THE GOVERNING EQUATIONS 


The equations of motion for isotropic media, in terms of the displacement u are [I]: 

(X. + |i) Uj ii + \lUjjj = p ti t (1) 
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or, in vector notation [1]: 


TV 

(\ + li)Vdiv(U) + iiV 2 U = p~f (2) 

d/ 

where X. and \l are Lame constants, and p is the mass density. Examining solution 
components of the form 


U = Ae 


i (K X - ojf) 


( 3 ) 


where U = {u x , u 2i u 3 ) T y A = (a,, a i) T > & - (k ly k 2 , k^) 1 , and X = (x Xy x 2 , x^) T y Eq. 
(2) results in: 


- (A. + H) (A ■ K)K-\i\K\ 2 A = -pco 2 4 (4) 

where \K\ 2 = K K, and in a matrix form: 


( - pw 2 + p|xi 2 + (X + p)* 2 

(X + ji) k ] k 2 

(X + |l)*,* 3 


( X + p) k j &2 

— pCO + P]/l| + (X + fi)&2 

{ X + p) ^2^3 

A 2 

< (X + p)*,*, 

(X + p.) 

- pco 2 + p| K \ 2 + (X 4- p) /c 2 j 

< J 


Non-trivial solutions exist when the determinant vanishes, i.e., when 

[ (^. + 2|i)|£l 2 - pco 2 ] [ fJ. | JS1 2 — p co 2 ] 2 = 0 (6) 

This can be written as: 


[c-^lAT] 2 - co 2 ] [c 2 s \K\ 2 -w> 2 ] =0 


where 


X + 2p, 
4 p 


f 1 

ip 


( 7 ) 


In the bulk, real K values yield the two types of waves that possess the dispersion relations 


to 2 = c]\K\ 2 and ft) 2 = c]\K\ 2 (8) 

which determine the phase velocity of the propagating waves given by Eq. (3). Given co, the 
set of solutions corresponds to two spheres in K space given by the above Eq. (7). The 
propagation of energy is given by the group velocity: 

y *- v '-4 and (9) 

thus, waves can propagate in all directions with speeds c g and c g . Fig. 1 shows the 
intersection of the group velocity vector plane with the k x lc 2 plane. For a given direction, the 
wave propagates with group velocity described by the vector from the zero origin to the 
surface of the plotted sphere. The vector length is c gp for the p-waves, and c gs for the s- 
waves. In general, when a wave reaches a boundary, other waves characterized by Eq. (7) 
too can be generated. As shown above, the media can support propagation of waves in any 
direction, thus the waves can propagate back into the interior. If this interface is an artificial 
boundary only, such waves are not acceptable. The purpose of a non-reflecting boundary 
region is to disallow such waves. 
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Fig. 1 



Group velocity vector in an isotropic solid as a function of K . 


EXISTING METHODS 


The treatment of artificial boundary can be divided into several main approaches; 

Viscous Boundary Conditions — Viscous boundary conditions have been proposed by 
Lysmer and Kuhlemeyer [2]. They have the form: 

du\ du 2 

a o c pJ t = a n and = a .2 (10) 

where p is the mass density, c p is the P-wave speed, c, is the S-wave speed, and w,, u 2 are 
the normal and the tangential directions to the boundary. The quantities a and b are 
dimensionless parameters, fitted to minimize the reflected energy for an incident wave at a 
particular chosen angle of incidence. This boundary condition was found to yield large 
spurious reflections at sharp angle of incidence. Some improvement of this boundary 
condition has been done [3]. This method is used in some NDE applications (see for 
example [8]). 

Non-reflecting Boundary Conditions — Absorbing boundary conditions for elastic 
waves have been developed by Clayton and Engquist [4], analyzed by Engquist and Majda 
[5], and improved by Higdon [7] who used the general absorbing boundary condition: 

i= '- 2 - 3 <"> 

where the number of terms is controlled by m. The appropriate choice of the P terms leads 
to absorption of different P- waves or S-waves at different angles of incident. The more 
terms are taken in Eq. (11), the less reflections are generated. Emerman and Stephen [6] 
showed that the original boundary conditions [4] are unstable for c s /c p < 0.46. 

Attenuating Buffer Zone — Cerjan et al [9], Sochacki et al [10], and Hanson & 
Petschek [11] presented buffer-domain techniques where the amplitudes of the 
displacement are gradually reduced when approaching the boundary. 

Non-local boundary conditions — have been studied by several authors and are 
recommended by Givoli [12]. These methods, however, are hard to implement and in 
general, low in computational efficiency. 
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NON-REFLECTING BUFFER 


The method presented in this paper uses a buffer zone as in the attenuating buffer- 
zone method. However, instead of attenuating the waves, we modify the equations of 
motion in the buffer in such a way that their direction of propagation is biased toward the 
outer boundary until no backward propagation is possible. 

The construction of the modified equations can be explained in terms of the dispersion 
relation. Consider a modified form of each of the factors of the determinant in Eq. (7) 

[c\K\ 2 - ( co - V 0 • AT) 2 ] = 0 (12) 


where V 0 = ( v,, v 2 , v ,) T . It yields the dispersion relation: 

(0 = c\K\ + V 0 K 


(13) 


and the group velocity: 


v*® = 4 + y 0 d4) 

Thus, the group velocity of the modified system can be controlled by an appropriate choice 
of Y 0 . Since K/\K\ is a unit vector, any Y 0 that satisfies 

V 0 m>max(c) (15) 

m being a direction vector, does not allow propagation of waves that have a component in 
the direction -m [Fig. 2], and the group velocity c p in the material is an appropriate choice 
for max(c). 

The change in the elasticity equation of motion that produces the above effect on 
velocity of the waves is given by: 


(X+*l) V div(U) + ]lV 2 U = p(| f + y 0 - V)V (16) 

When V Q is zero, Eq. (16) reverts to the physical Eq. (2). We see that we can use 
Y 0 = 0 in the physical domain of interest, and a non zero value of Y 0 satisfying Eq. (15), at 
the outer boundary of the buffer region, where m = n, n being the outward normal at that 



Fig. 2 Group velocity in the non-reflecting buffer. 
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boundary. As a discontinuity in V 0 in the buffer zone will produce reflection by itself, V 0 
must change smoothly from zero to its maximum value. 

The analytical study of the change in V 0 and the amount of reflected energy as a 
function of Y 0 will be presented elsewhere. Certainly, the larger the buffer domain is, the 
less reflection it will produce by an appropriate smoothly varying V 0 . 


SURFACE WAVES 


In order to analyze the behavior of the surface waves in the buffer zone, we consider 
waves of the form of Eq. (3) with K = (k ]y 0, ifc 3 ) , where k u k 3 are real, i.e.. 


-itt)/ ik.x 

e e e 


(17) 


which are components of surface waves in the ±x { direction, propagating on the z = 0 
surface. For simplicity of the analysis we consider the outward normal to the nonreflecting 
boundary zone to be n - (- 1 , 0, 0) , and V 0 = (-v, 0, 0) , i.e., V 0 - K = -vk ] . 
Substituting these K and J/ 0 into Eq. (12), the resulting dispersion relation is 


c (k]~ +vky = 0 


for 


c c 


(18) 


Rewriting it as a quadratic equation in &,/ 0) we get 


2 i k] k y 

( c - v ) ( — ) -2v— - 
03 (0 


2 k, 

1+^ (-) 
0) 


= 0 


(19) 


The product and the sum of the two roots of Eq. (19) £| (l) /0) and (fc, U) / to) are: 


, O) , (2) 1 + C (— ) 

kt k i CO 


CO CO v 2 - c 2 



CO CO 



( 20 ) 


The requirement v > c from Eq. (15) implies that the real parts of the two roots are negative 
and the waves propagate in the direction of n, i.e., waves do not propagate backwards. It can 
be shown that V 0 • n must be non-negative at any boundary point, including any physical 
boundary, otherwise, surface waves will be amplified in an uncontrolled way. 


DISCRETIZATION 

The construction of the buffer region equations as described above is in the 
continuous level, independent of the discretization. For the numerical solution, a staggered 
grid discretization in space was used [ 131. A central second order scheme was used for all 
the spatial terms of Eq. (16) except for the V 0 • V term. This term, which exists only in the 
buffer zone, requires careful discretization. Considering the resulting dispersion relation for 
the discrete system, we can see that a central differencing, for example, is inappropriate. 
This can be seen by analyzing the dispersion relation. The discrete version of Eq. (3) is: 


r7 . HQ X/h- <of) 

U = Ae 


( 21 ) 
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h being the mesh size. The dispersion relation, using a central scheme for all the terms, 
including F 0 • V, derived in an analogous way to the calculation of Eq. (6), is 


0 9 0 

4(X. + 2p) (sin 2 -^ + sin 2 y + sin 2 y )-p (ca-v^inG^VjsinOj-VjsinGj) 2 
- 0 0 0 2 ( 22 ) 
• 4p(sin 2 y + sin 2 y + sin 2 y)-p (co-v^inG,— v 2 sin0 2 -v 3 sin0j) 2 l =0 


It can be seen that for real ( 0 , if (e,,0 2 ,0 3 ) is a solution, then (±0,, ±0 2 , ±0 3 ) is also a 
solution. A change of sign in 0 . corresponds to a wave propagating in the i direction 
opposite to the original wave. Thus, although the modified differential equation produces no 
reflections, some of its discretizations might do so. 

This is resolved by using an upwind discretization for V ■ V. It changes the 
(-v^inG^VjsinGj-VjSinGj) term in Eq. (22) to 


-v, ( 1 - e“ 9 ') -v 2 ( 1 - i*') -v 3 ( 1 - e"*’) (23) 


for v, > 0, v 2 > 0, v 3 > 0. It can be shown that the resulting dispersion relation does not allow 
reflected waves when V 0 ■ o is larger than the maximal group velocity in the discrete system, 
which for the above discretization is c„ . 


APPLICATION 


The method was tested on a body with rectangular geometry. An artificial buffer 
domain that was 32 grid points wide was positioned at the left surface, while the other 
surfaces where physical ones, with stress-free boundary conditions. The result in Fig. 3 
shows the horizontal displacement for a body under plane strain conditions due to a 
localized dynamic force acting on the top surface. We can see the bulk p and s waves, as 
well as the surface waves propagate as a function of time. Mode conversion at the boundary 
is clearly visible. As the waves enter the buffer region, they speed up and attenuate due to 
the upwind discretization for the added V ■ V term. No reflections occur at that artificial 
boundary, while at the same time, the waves that reach the physical boundary on the right 
are reflected as expected. Note that even a comer between an artificial boundary and a 
physical one does not produced reflected waves. 


CONCLUSIONS 


The proposed method for constructing a non-reflective region to serve as an artificial 
boundary is easy to apply, and depends only on the maximum group velocity of the 
discretized system. It is very effective for practically eliminating reflection of p-waves and 
s- waves as well as surface waves. The method is independent of the discretization and can 
be used with high order discretization as well. 
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Dynamic load 









Fig. 3 Horizontal displacement due to a localized dynamic load acting on the top 
surface. The left boundary includes a non-re flective region. The other boundaries are 
stress-free. 
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